Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 18 February 2013 (MN style file v2.2) 



The Neutral Hydrogen Content of Galaxies in 
Cosmological Hydrodynamic Simulations 

Romeel Dave^'^, Neal Katz^, Benjamin D. Oppenheimer^, Juna A. Kollmeier^, 
"cn David H. Weinberg^ 



O 

CD 



o 

u 

o 



cn 
m 

O 

m 



University of the Western Cape, 7535 Bellville, Cape Town, South Africa 
Astronomy Department, University of Arizona, Tucson, AZ 85721, USA 
Astronomy Department, University of Massachusetts, Amherst, MA 01003, USA 
Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, Netherlands 
Observatories of the Carnegie Institute of Washington, Pasadena, CA 91101, USA 
Astronomy Department and CCAPP, Ohio State University, Columbus, OH 43210, USA 



18 February 2013 



ABSTRACT 

We examine the global H i properties of galaxies in quarter-billion particle cosmologi- 
cal hydrodynamic simulations, focusing on how our main adjustable physical process, 
galactic outflows, impacts H i content. In addition to the three outflow models consid- 
ered in our earlier papers, we present a new one (ezw) motivated by high resolution 
interstellar medium simulations, in which the scalings of wind speeds and mass loading 
factors follow those expected for momentum-driven outflows for larger galaxies, and 
energy-driven outflows for dwarf galaxies (cr < 75 kms~^). To obtain predicted Hi 
masses, we employ a simple but effective local correction for particle self-shielding, as 
well as an observationally-constrained transition from neutral to molecular hydrogen. 
We find that our ezw model produces an H i mass function whose shape agrees well 
with observations from the ALFALFA survey, having a low mass end slope of —1.3, 
while other models agree less well. Outflows critically govern the H i content in low- 
mass galaxies, with higher mass loading factors yielding higher Hi fractions. Satellite 
galaxies have a bimodal distribution in Hi fraction versus halo mass, with lower mass 
satellites and/or satellites in larger halos more often being devoid of Hi. At a given 
stellar mass, H i content correlates with star formation rate and inversely correlates 
with metallicity, as expected if driven by stochasticity in the accretion rate. At higher 
redshifts, massive H i galaxies become less frequent and the H i mass function becomes 
significantly steeper. The global cosmic Hi density conspires to remain fairly constant 
from z ~ 5 ^ 0, but the relative contribution from smaller galaxies increases substan- 
tially with redshift. Overall, Hi in galaxies reflects a transient reservoir of fuel for star 
formation, and hence provides a crucial glimpse into the inflow and outflow processes 
that govern galaxy evolution. 



1 INTRODUCTION 



The current paradigm for galaxy evolution rests on the tenet 
that gas flows into and out of gala^cies are primarily respon- 
sible for go v ernin g galaxy growth (e.g. [iCcrcs ct al . 2005, ; 
iDekel et al.l l2009l : iDave. Finlator. fc Qppenheimeil |2012| . 
and references therein). Such accretion and outflow pro- 
cesses are often difficult to detect directly owing to the ten- 
uous and multi-phase nature of the gas, but they are ex- 
pected to leave clear imprints on the stellar and gaseous 
content of galaxies. A key region where such baryon cycling 
processes might be probed is in the outskirts of star-forming 
galaxies, where reservoirs of neutral hydrogen (Hi) are ex- 
pected to hold the raw material that ultimately fuels star 
formation. This reservoir is thought to be continually re- 
plenished by gravitationally-driven accretion from the highly 



ionised intergalactic medium (IGM), possibly augmented by 
re-accretion of ejected gas. 

Observations of the neutral hydrogen content of galax- 
ies using 21cm fine structure line emission are progressing 
rapidly, and promise to improve further. The Hi Parkes All- 
Sky Survey (HIPASS; Me yer et al. 2004) provided a com- 
prehensive deep census of the Hi content of low-redshift 
galaxies. This was f ollowed by the Arecibo Fast Legacy Alfa 
Survey (ALFALFA: iGiovanefli et al.ll2005l ). which surveyed 
7000 deg^ for Hi 21cm emitting objec ts. More re cently, the 
GALEX Arecibo SDSS Survey (GASS: ICatinella e t al. 201^) 
has obtained high-fldelity Hi data on optically-selected 
galaxies, to provide a more Hi-unbiased census. Construc- 
tion is underway of the MeerKAT telescope, which will pro- 
vide the ability to detect H I in galaxies out to unprecedented 
redshifts (e.g. the Looking At the Distan t Universe with 
the MeerKAT Array - LADUMA - Survey: iHolwerda et al] 
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[ioil), and the Australian SKA Pathfinder (ASKAP), which 
will analogously perform the Wallaby Survey. Ultimately, 
the Square Kilometer Array (SKA) will enable studies of 
Hi both nearby and out to higher redshifts that far exceed 
today's capabilities. Complementary information is provided 
by studies of damped Lya absorbers (DLAs) that trace 
Hi in and a round galaxies in t he spectra of background 
quasars (e.g. iBattisti et"alll2012l ') . 



It is therefore timely and crucial to develop a theoret- 
ical framework for understanding the physics that governs 
the H I content of galaxies and its evolution. Because galac- 
tic H I reservoirs are thought to represent a transient phase 
of baryons (jProchaska & Wolfe 2 009.) as they pas s from the 
diffuse intergalactic medium (e.g. lOave et al.f2010l ) to higher 
density molecular gas, which eventually converts into stars, 
it is essential that models for cosmic H l be dynamical in na- 
ture. Cosmological hydrodynamic simulations provide one 
such type of dynamical model for this purpose, allowing one 
to directly track the physical state of gas as it flows through 
the Hi reservoir. Particularly, the inclusion of galactic out- 
flows in simulations provides some interesting new twists 
on the accretion paradigm, as enriched outflows that return 
back into galaxies may p rovide a significant sourc e of ac- 
cretion ( "wind recycling" ; lOppenheimer et al.ll2010l ). These 
inflow, outflow, and recycling processes are expected to man- 
ifest themselves in the Hi content of galaxies. 



2 METHODS 



2.1 Simulations 



Our simulations are evolved with an extended version of 
the Gadget- 2 N-body -|- Sm oothed Particle Hydrodynamic 



In this work, we examine the Hi content of galaxies 
in cosmological hydrodynamic simulations with a variety of 
outflow models . Thi s work follows on our early study in 
[Popping et al.l (|2009h . but it uses a significantly improved 
simulation suite and focuses on the physical processes that 
govern the Hi content in addition to making predictions 
for Hi observables and t heir evolution. This work is also 
comparable to that of Duffv et al.l (|2012l ). who used simu- 
lations to examine both the atomic and molecular content 
of simulated galaxies. Compared to that work, our model 
explores a different range of feedback parameters (includ- 
ing ones that match Hi observations substantially better), 
probes down to significantly smaller systems (albeit in a 
smaller volume), and uses a slightly different approach to 
compute the Hi content of galaxies. Nonetheless, for over- 
lapping mode ls and mass r anges , our results generally agree 
with those of iDuffv et al] l|2012l) . 



The paper is organised as follows. In Sj2] we introduce 
our simulations, and describe our method for calculating the 
atomic and molecular gas content of our simulated galaxies. 
In 33] we examine the H I mass function and its evolution, 
and we examine the complementary constraint of H I rich- 
ness versus stellar mass in SjH In [J5]we study how H I content 
is impacted by environment, and in !j6]we show how it re- 
lates to galaxy metallicity and star formation rate. In iJ7|we 
make predictions for the evolution of Q.m, and compare to 
observations out to z ~ 4. In ij8]we examine key Hi prop- 
erties in a simulation with different numerical resolution to 
test for convergence. Finally, we summarise and discuss our 
results in f(9l 



dy - 

fSP H) code (ISpringe]||2005l ). We assume a ACDM cosmol- 



ogy (|Hinshaw et al.ll2009l ): f^M = 0.28, = 0.72, h = 
Ho/{lQQ kms~^ Mpc~^) = 0.7, a primordial power spec- 
trum index n = 0.96, an amplitude of the mass fluctua- 
tions scaled to (Tg = 0.82, and fib = 0.046. We call this 
cosmology our r-series, where our general naming conven- 
tion is T[boxsize]n[particles/side][wind model]. Our primary 
simulations use a cubic volume of 32/i~^Mpc on a side with 
512"^ dark matter and 512^ gas particles, and a softening 
length of e = 1.25^~^kpc (comoving, Plummer equivalent). 
The gas particle mass is 4.5 x lO^i\f0, and the dark mat- 
ter particles masses are approximately five times larger. We 
can thus reliably resolve galaxies down to stellar masses of 
Af*,iim = 1.4 X 10* (see discussion below). 

Our version of Gadget-2 includes cooling p rocesses us- 
ing th e primordial abundances as described by iKatz et al.l 
(| 19961 ). with additional cooling from metal lines a s sumin g 
photo- ionisation equilibrium from Wi ersma et al.l (|2009l ). 
St ar formation is modelled using a subgrid recipe introduced 
bv lSpringel fc HernauistI diool ) where a gas particle above 
a density threshold of nn = 0.13 cm~^ is modelled as a 
fraction of cold clouds embedded in a warm ionised medium 
following McKcc & Ostrikcr (1977). Star formation (SF) fol- 
lows a Schmidt law (|Schmidtj|l959l ') with the SF timescale 
scaled to matc h the z = Ke nnicutt relation jKennicuttl 
1998). We use a lChabriej (|2003l ') initial mass function (IMF) 
throughout. We account for metal enrichment from Type II 
supernovae (SNe), Type la SNe, and AGB stars, and we 
tr ack four elements (C,0,Si,Fe ) individually, as described 
bv lOppenheimer" fc Dav j (|2008l ). 

Galactic outflows are implemented using a Monte Carlo 
approach analogous to star formation. Outflows are directly 
tied to the SFR, using the relation M^i^d ~ JyxSFR, where 
rj is the outflow mass loading factor. The probability for a 
gas particle to spawn a star particle is calculated from the 
subgrid model described above, and the probability to be 
launched in a wind is rj times the star formation probabil- 
ity. If the particle is selected to be launched, it is given a 
velocity boost of Vw in the direction of vxa, where v and 
a are the particle's instantaneous velocity and acceleration, 
respectively. Once a gas particle is launched, its hydrody- 
namic (not gravitational) forces are turned off until either 
1.95 X 10^°/(ww(kms-^)) years have passed or, as more of- 
ten occurs, the gas particle has reached a density that is 10% 
of the SF critical density (i.e., 0.013 cm~^). This attempts to 
mock up chimneys generated by outflows that would allow a 
relatively unfettered escape from the galactic ISM, a process 
not properly captured by the spherically-averaging SPH al- 
gorithm at ~kpc resolution. It also yields results that are 
less sensitive to numerical resolution jSpringel fc HernauistI 
2003b) than models that do not turn off hydrodynamic 
for ces. For a further discussion of h ydrodynamic decoupling, 
see lDalla Vecchia fc Schav^ (|2008l ). 

Choices of the parameters 77 and define the "wind 
model". For this paper we make use of the following four 
wind models: 

(i) No winds (nw), where we do not include outflows 
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(i.e. ?7 = 0). This model fails to match a wide ra r iRe of 
observables (e.g. Dave. Qppenheimer. fc Finlatorl I2OIII : 



iDave. Finlator. fc OppenheimeJ 120111 1. but is included to 
establish a baseline for the overall impact of winds. 

(ii) Constant winds (cw), where -q = 2 and 
Ww = 680 kms^^ for all galaxie s. This mo d el is simi- 
lar to the wind model employed bv lDuffv et al] (|2012h . 

(iii) Momentum-conserving winds (vzw), 
where the wind speed and the mass loading fac- 
tor depend on the galaxy velocity dispersion 
o" (^M urrav. Quataert. fc Thom pson 2005j), scaling as 
estimated (see lOppenheimer fc DavQi2008l ). 



= 3(T-^/ /i - 1, 



V 



(1) 
(2) 



where /l ~ [1.05, 2] is the luminosity factor in units of 
the galactic Eddington luminosity (i.e. the critical lumi- 
nosity necessary to expel gas from the galaxy potential), 
and (To = 150 km/s is the normalisation of the mass load- 
ing fact or. Choices for the former are t aken from obser- 
vations l|Rupke. Veilleux fc Sanders! |2005| ). while the latter 
is bro adly constrained to match h igh-redshift IGM enrich- 
ment (|Oppenheimer fc Dave'2008'V Galaxies are identified 
"on-the-fly" during the simulation using a Friends-of- Friends 
(FOF) algorithm applied to the gas, star, and dark matter 
particles. We choose a smaller than usual linking length to 
pick out just the galaxies. This linking length evolves with 
redshift and in terms of the mean interparticle separation of 
all particles is 



0.06 



Ho 



(3) 



We estimate the velocity dispersions necessary for our wind 
scaling laws using the total FOF galaxy mass MgalU using 
the relation 



200 



gal 



5 X 10i2/i-iA'/<7, Ho 



km/s. 



(4) 



Which we have empirically determined gives an accurate 
measure of the velocity dispersion in our simulations. 

(iv) Hybrid energy /momentum-driven winds 
(ezw), which employs the vzw scalings for galaxies with 
o" > 75 kms~^, then switches over to a steeper depen- 
dence of 77 oc a"^ in cr < 75 kms^^ systems. The wind 
speed still scales proportionally to a as in the vzw model. 
This model roughly captures the behaviour in recent ana- 
lytic and hyd rodynamic models of outflows from inte rstel- 
lar media bv iMurrav. Quataert. fc ThompsonI (|2010l ) and 
iHopkins. Quataert. fc Murravl ( 20121) . respectively. The ba- 
sic idea is that in dwarf galaxies the energy from super- 
novae plays a dominant role in driving outflows, while in 
larger systems the momentum flux from young stars and/or 
supernovae is the dominant driver. As a result, the out- 
flow scalings switch from momentum-driven at high masses 
to energy-driven at low masses. We make this transition 
abruptly at ct = 75 kms^^ guided by the analytic models 



^ We hav e actually been using this alg orithm, which slightly dif- 
fers from llOppenheimer in all of our recent simu- 
lations starting with JOppenheimer "et aLlbOlOh . 



of iMurrav, Quataert. fc Thom pson' ("2005', '20ld ) and high- 
resolu tion ISM simulations of iHopkins. Quataert. fc Murravl 
l|2012f ). although it should perhaps be more gradual (note 
that ?7 itself is continuous across this boundary). In any case, 
this model captures the gist of the most up-to-date small- 
scale outflow models. As we will see, the ezw outflow model 
fares somewhat better compared to observations than vzw, 
which in turn fares much better than cw or nw. 

Additionally, in our fiducial ezw simulation, we employ 
a heuristic prescription to quench star formation in massive 
galaxies. This is not a physical model, but simply a tuned 
parametrisation to limit star formation in massive systems 
to reproduce the observed exponential high-mass cutoff in 
the stellar mass function and to make the simulation faster 
to run. In this prescription, we quench star formation in an 
entire galaxy according to a probability, Pq, given by the 
equation 



2 L log(crQspr) 



(5) 



where the median a at which a galaxy has a 50% chance 
of being quenched is aqmcd ~ 110 kms~^, corresponding 
to A/halo = 10^^'^ Mq at z — 0, and the spread in a is 
CQspr ~ 32 kms~^. We also require ct > 75 kms^^ to fully 
suppress the already low probability of lower mass galax- 
ies being quenched. We have found that these parameter 
choices nicely reproduces the high-mass end of the stellar 
mass function, as we will show in i]2.3l 

Every time we identify galaxies using our FQF group 
finder, which we do to calculate a for the vzw and ezw wind 
models and occurs approximately every 10 Myr, we compute 
this quenching probability for each identified galaxy. If the 
galaxy is "quenched", then each time a gas particle would 
have formed a star over the time interval until the next time 
we identify galaxies using our FOF group finder, it is instead 
heated to 50 times the virial temperature, Tvir, where 

, 2 



Tvir = 5 X 10" 



K. 



(6) 



.200 kms-i. 

The motivation for heating the gas to such extreme temper- 
atures is primarily to prevent it from re-accreting at later 
times and thus requiring multiple ejections; the total en- 
ergy input is thus less for higher quenching temperatures. 
The median temperature to which the ISM is heated by 
quenching is 10*'^ K and arises from a median halo mass 
of 10^^'^ M0. The energy input from quenching averages to 
10*'''^ergs~^Mpc~"^ between z = 0.75 — 2.5, corresponding to 
the peak of AGN activity, and the integrated energy input 
until z = is 8x 10^* ergs/g, which equals 9x 10~^ of the rest 
mass energy of all baryons. Considering that 6.5% of baryons 
are in stars in this simulation at z = 0, and using the as- 
sumption that supermassive black holes (SMBHs) have 10~^ 
of the mass in stars, the quenching energy corresponds to 
1.4% of the rest-mass energy of SMBHs, which is compara- 
ble to the energy imparted from AGN feedback in cosmologi- 
cal simulations that self-consistently track black hole growth 
and feedback (Pi Matteo et al. 2005: B ooth fc Sc havc 200^). 

In this paper we are mostly concerned with galaxies 
around L* and below, which are mostly unquenched. We 
emphasise that this quenching prescription is not intended 
to be a realistic physical model, and is actually in large part a 
computational convenience, as removing gas from the largest 
galaxies substantially speeds up the code at low redshifts 
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and allows us to run our large simulations to z = within 
practical time frames. We refer the reader to iDuffv et alj 
l|2012l ) for a discussion of Hi in larger- volume simulations 
that attempt to quench massive galaxies based on directly 
tracking black hole growth and feedback. 

Finally, our no-wind simulation employs a volume of 
16/i~^Mpc with 2 X 256'^ particles. This results in the same 
mass and spatial resolution as the wind simulations, but in 
a volume that is 8 times smaller. As we mostly use the no- 
wind simulation to qualitatively show the impact of winds 
at small masses, the reduced volume for simulation will not 
significantly affect our conclusions. 



2.2 Computing the Hi content 

We use SKIE0 (Spline Kernel Interpolative Denmax) to 
identi fy galaxies as boun d groups of star-forming g as and 
stars teeres et al.) l2005l : lOppenheimer et al.l I2OI0I ). Our 

gala xy stellar mass lim it is set to be ^ 64 star parti- 
cles (|Finlator et al.|[2006l ). resulting in a minimum resolved 
mass of M.^min ~ 1.4 X 10* Mq. We will only consider galax- 
ies with stellar masses M« ^ M,,min in our analysis, regard- 
less of their Hi content. Our resolution convergence tests in 
38]show that the H I properties are reasonably well converged 
even at this stellar mass threshold. 

We identify dark matter halos via a spherical overden- 
sity algorithm (Kcrcs ct al. 2005) out to a virial overdensity 
given m fPave et al.l (|201d . see their eq. 1), centred on each 
galaxy. We separate galaxies into central and satellite galax- 
ies by associating each galaxy with a halo; if a galaxy's cen- 
tre lies within the virial radius of a larger (by stellar mass) 
galaxy, we consider it to be a satellite of that galaxy, and 
the halos of those two galaxies are merged. Note that for 
galajcies near the edge of a larger halo, this can result in 
halo that has "bumps" along its (mostly spherical) surface. 

To compute the Hi content of galaxies, we need to iso- 
late the gas that is in the (predominantly) neutral phase. 
This involves defining two "boundaries": the division be- 
tween gas that is exposed to the full metagalactic ionising 
flux and gas that is self-shielded, and the division between 
atomic and molecular gas. 

To calculate the self-shielding, we first compute the 
neutral hydrogen component of each gas particle under 
the assumption tha t it is not self-shielded. We follow 
[Popping et al.l l|2009t ) who employed a simple hydrogen ioni- 
sation balance calculation that yielded the following formula 
to determine the neutral fraction: 



2C + 1 - V(2C + 1)2 



4C2 



with 



C - 



2C 



np{T) 



(7) 



(8) 



where n is the Hydrogen number density, T is the gas 
temperature, Fhi is the Hi photo-ionisation rate, and the 
recombination rate coefficient /? is well-fit by the func- 
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Figure 1. Gas particle properties as a function of number density 
riH in our r48n384vzw simulation. Left panels sliow z = 0, right 
panels show 2 = 3. Top: Hi fraction (blue) and H2 fraction (red) 
based on the calculation described in the text. We show a random 
sampling of 1% of all the particles. Middle: Total mass density 
Q per unit lognn in Hi (blue) and H2 (red). Bottom: Pressure- 
density relation. Green points show unshielded particles (either 
below the pressure threshold or above the density threshold), blue 
points show Hl-dominated gas, and red points show molecular- 
dominated gas. 



tion ijVerner fc Ferlandlll996l l 

/?(r) = a[y(r/r„)(i + ^iT/To)f-\i + \/(t/Ti))'+'] 

, (9) 

For H I, the fitting parameters are a — 7.982 x 10 cm^s ^ , 
b = 0.7480, To = 3.148 K, and Ti = 7 036 x lO'' K. We take 
Fhi from the lHaardt fc Madau I (|200lh ionising background, 
whose amplitude is adjusted to match the observed mean 
flux decrement i n the Lya forest in these simulations (see 
iDave et al.ll2O10l ): we describe this further below. 

With the optically-thin neutral fraction computed for 
each gas particle, we employ a simple particle-by-particle 
post-processing correction for self-shielding of the Hi. We 
resort to this approach because it is computationally pro- 
hibitive to do the full radiative line transfer on these sim- 
ulations; instead, we will c alibrate our approach to radia- 
tive tr ansfer simulations by iFaucher-Giguere. Keres. fc Mai 
ll2010ll . 

We begin by assuming that each particle has a density 
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profile given by tiie SPH spline kernel W{r) (see ISpringej 
[2005, for definition). We compute the radial column density 
profile as follows: 



iVHi(r) = 



0.76/HIP9 



W(r')dr', 



(10) 



where pg is the SPH density of the gas particle, nip is the 
proton mass, and h is the particle's smoothing length. We 
then determine the radius R where Nhi{R) ~ A'ni.iim, with 
A'ni.iim = lO^^'^cm"^ being where the particle becomes op- 
tically thick (r — 1) to continuum photons at the Lyman 
limit. If no such radius R exists, then the particle is opti- 
cally thin. If R exists, then we compute the unshielded and 
shielded mass fraction using 47r J W{r)r^dr integrated from 
h ^ R and 7? — >■ 0, respectively. The neutral fraction of 
the particle is then the mass-weighted mean neutral frac- 
tion, assuming that the portion from _R — >■ is 90% neutral 
(since some ionised gas is seen to exist even within dense re- 
gions). By this procedure, particles in dense regions "auto- 
shield" themselves from the ambient flux, though the effects 
of shielding from nearby gas are neglected. 

This computation relies on knowing the metagalactic 
photo- ionising flux, since this determines the neutral frac- 
ti on in the optically - thin regime. Our simulations assume 
a lHaardt fc Madau I (|200ll ) ionising background, but more 
detailed constraints can be placed on the amplitude of the 
Hi photo-ionisation rate by using the measured mean flux 
decrement in the Lyman alpha forest. In particular, at each 
redshift for each simulation, we determine a "flux factor", 
which is the value by wh i ch we must multiply the strength 
of the lHaardt fc Madau I (|200lD background to achieve con- 
sistency with the observed flux decrement. For the observed 
flux decrement , we em ploy at z < 2 the determination by 
iKirkman et al.l ll2007l'l. while for z ^ 2 we take the val- 



ues from Becker et al. (|2012h . As described more fully in 



iDave et all l|2010l ). we extract 1000 spectra and iteratively 
adjust the flux factor until the mean flux is within 1% of the 
observed value. The resulting flux factors for our ezw sim- 
ulation at z = 0,1,2,3,4,5 are 1.72, 2.03, 1.66, 0.92, 1.13, 
and 1.70, respectively. For our vzw and cw simulations, the 
z = flux factors are 1.63 and 1.41 respectively (we do not 
consider these models at higher redshifts). 

Next, we separate the molecular component from the 
atomic Hi. To do so, we employ the observed ISM pres- 
sure relation from The H I Nearby Galaxy Survey (THINGS; 
iLerov et al ]|2008l ): speciflcally, we use their "combined spiral 



subsample" 
gas as 



flt, which gives the ratio of molecular to atomic 



Rmol = (P/PoT, (11) 

where Po = 1-7 x 10"' cm~^K and a = 0.8. We compute the 
gas pressure based on the density and two-phase medium 
temperatur e of each star-forming gas p article (our prescrip- 
tion follows ISpringel fc Hernauistll2003l '). Note that we only 
compute molecular fractions for star-forming gas particles, 
which in our simulations is gas with nu > 0.13 cm~^. Gas 
that is not star-forming is assumed to have zero molec- 
ular content, and in any case the formula above would 
yield a very small mol ecular fraction. Thi s prescription fol- 
lows that emp l oyed by iDuffv et al.l (|2012l '). but differs from 
iPopping et al] (|2009l ) who used a fixed pressure threshold of 
810 cm~^K to separate atomic from molecular gas. 



Figure[T]illustrates the resulting gas particle atomic and 
molecular fractions. The top left panel shows the H I and H2 
mass fractions in particles from our momentum-driven wind 
scalings (vzw) simulations at 2: = 0. In the optically thin 
regime (mostly not depicted here; see lDave et al.|[2010l . for 
a full phase space diagram), /hi scales linearly with nn- 
Then there is a sharp upturn at nu ~ 2 x 10~'^cm~^, when 
auto-shielding becomes important. 

At high redshifts, the stronger ionising background 
causes auto-shielding to set in at a higher density. A sim- 
ple scaling shows that A^Hi.Um oc fmnnl, where the length 
/ oc n^^^^ for a spherical cloud (or SPH particle), and 
/hi oc nu/Tm; hence A'ni.iim oc /Tm- Therefore, a fac- 
tor of 10 increase in Phi, which occurs between z = — 3, 

will lead to a factor of 4 increase in the nn where self- 

3/5 

shielding becomes important, because tih oc Pjjj for con- 
stant A'^Hi.iim. Figure [1] (top right panel) shows that auto- 
shielding at z = 3 is effective above uh ~ lO^^cm"'^ as 
expected for the stronger z — 3 background UV field, which 
also matches the z = 3 radia t ive tr ansfer simulations of 
Fauchcr- Gisuere^, Keres. fc Mai (|2010l , see their Figure 3). 
Indeed the overall shape of fm{nH) is actually quite sim- 
ilar to theirs, although there is less scatter at a given uh in 
our prescription owing to the fact that we do not consider 
shielding from neighbouring particles. This shows that our 
physically-motivated choice of auto-shielding with A^'Hi.iim 
at the Lyman limit is a reasonable approximation to much 
more detailed radiative transfer models. 

Moving to higher densities, eventually we reach the 
star-forming density threshold, above which the molecular 
fraction becomes nonzero. At hh ^ 0.5 cm~^, correspond- 
ing to a pressure close to Po, the gas becomes molecular- 
dominated. The pressure relation is shown in the bottom 
panels of Figure [TJ with the particles colour-coded by their 
dominant phase of hydrogen. The THINGS pressure thresh- 
old is_jnu£h_Jngher than the P/k — 810 cm~''K assumed 
in [Popping et al.l l|2009t ) and, moreover, our prescription 
produces a more gradual transition between atomic and 
molecular-dominated gas. A change in slope at nn = 0.13 
cm""' occurs because this is the density above which we al- 
low star to form and the gas particles become two-phase. 

In this paper we are mostly interested in the regime 



10" 



< Uh 0.5 cm where the majority of cosmic H I 



resides (at z = 0), above which gas becomes mostly molec- 
ular and below which it is mostly ionised. This is shown 
in the middle panels of Figure [TJ where we plot the total 
mass density f2 in Hi (blue) and H2 (red) per unit inter- 
val of lognn; the peak is around uh ~ 10^^ cm~^. This 
Hl-dominant density regime is reasonably well resolved in 
our simulations, and hence the predicted H I content is ex- 
pected to be robust, despite it being a transitory phase. In 
contrast, the transition from molecular gas to stars typically 
occurs at densities well above what we can resolve directly, 
and hence the molecular content will not be as robustly pre- 
dicted; for this reason comparisons to the molecular con- 
tent of galaxies may be less robust (thoug h still interesting; 
we refer the reader to iDuffv et al.l [2OI2I for such predic- 
tions). Nonetheless, our predictions for the evolution from 
z ~ 2 — >■ of the total star-forming gas (most of which is 
typically molecular) in our momentum-driven wind (vzw) 
simulation are in good agreement with observations from 



© 0000 RAS, MNRAS 000, 000-000 



6 Dave et al. 




CO 

O 

Oh . 



o ■ 



9 10 11 

log M. (MJ 

Figure 2. Galaxy stellar mass functions at z = in our four 
wind models: ezw (green), vzw (cyan), cw (red), and nw (blue). 
The vertical dotted line is the stellar mass resolution limit. 
The thick solid line shows a fit to the observed GSMF from 
iBaldrv. G lazebroo k. fc Drivej ||2008|V The ezw outflow model, 
which includes a heuristic prescription for quenching star for- 
mation at high masses, matches the observed GSMF to within 
uncertainties. 



the Plateau de Bure High- 2: Blue Sequence Survey (PHIBSS; 
iTacconi et al. I [20121 . see their Figure 13); other wind models 
fare less well. Finally, we note that the stellar mass growth 
rate is fairly robust because it is not a transitory phase but 
an end state of accre ted gas, an d is typicall y limited by 
the gas supply rate (e.g Finlator fc Dave 2008; Bouche et al.l 
I2OI0I : iDave. Finlator. fc Oppenheimeill2012 ) . 

With each gas particle's Hi content determined, we 
must now associate the gas particles with galaxies. The in- 
formation from SKID is not sufficient, because SKID only in- 
cludes star-forming gas and stars in a galaxy, while a signifi- 
cant amount of the H I resides in an extended region around 
the galaxy, beyond the actively star-forming region. Thus 
to account for extended H I, we add up all the H I mass in a 
sphere around each galaxy that extends to the outermost ra- 
dius as defined by SKID, i.e. the radius of the farthest SKID 
particle associated with that galaxy. The outermost radius 
is typically many times the half-mass radius. While it may 
seem that this still may not fully account for an extended H I 
disk, in practise the low threshold density for star formation 
in our simulations means that this choice still encompasses 
the vast majority of the Hi. We tried extending this to 1.5 
times the outermost radius, and the total Hi mass in the 
volume increased by only 4%; furthermore, one will start to 
increasingly "double count" gas that may be between nearby 
galaxies as being part of both galaxies. 



The ezw model provides a strikingly good fit to the 
observed GSMF. At the massive end, this is a direct con- 
sequence of tuning the quenching prescription as described 
previously. In contrast, the low-mass end is unaffected by 
our quenching prescription, and instead reflects the effect of 
the hybrid energy-momentum driven winds. There is a clear 
upturn in the GSMF at M* < Iff'^Mp,, which is also seen 
in the data. As explained in lOppenheimer et al.l l|2010l ). in 
our simulations this arises because above this mass, wind re- 
cycling becomes increasingly important, and provides extra 
fuel to higher mass galaxies. Below this mass, the typical 
recycling time becomes longer than a Hubble time, and the 
slope begins to steepen towards the dark matter halo mass 
function's slope. 

The vzw model provides not quite as good a 
fit to the low mass end, similar to what wEis seen 
in a lower-resolution version of th e same model in 
iDave. Qppenheimer. fc Finlatoil (|201ll ). The differences at 
the massive end effectively show the impact of the quenching 
model, since our ezw simulation uses it while our vzw sim- 
ulation does not, while the two outfiow prescriptions them- 
selves are identical in this mass regime. Quenching has a 
substantial effect on the GSMF, but we will show that it 
has a minimal effect on the H I mass function. Meanwhile, 
the constant wind model produces a very steep low mass 
end slope that looks nothing like the data, while the no 
wind model overproduces stars at virtually all masses be- 
cause it strongly overcools baryons. These results again fol- 
low those obtained using lower-r esolution simulations of the 
sam e wind models presented in lOppenheimer et al.l (|2010l ) 



2.3 The Stellar Mass Function 

As an initial baseline statistic to compare our four wind 
models, we show in Figure [2] their galaxy stellar mass func- 
tions (G^MFs}_a^^_^_0;_VVe_£om to observations 
from lBaldrv. Glazebrook. fc Driver! l|2008l ) using SDSS. 



and iDave. Qppenheimer. fc Finlatorl ~ 20 111) , and more dis- 
cussions of the GSMF can be found there. 

This ezw simulation is the first hydrodynamic simula- 
tion that we are aware of to yield agreement with the ob- 
served GSMF to within statistical uncertainties over the en- 
tire mass range probed. This does not imply that this model 
is fully correct, as it could still be that the star formation 
histories in this model are incorrect, and other constraints 
may not be as well matched. As a case in point, we will 
show in ijS]that the mass-metallicity relation in this model 
looks too steep compared to observations, and that the vzw 
model provides a better match. Furthermore, the success of 
this model at high masses owes to including an ad hoc and 
highly tuned prescription for quenching star formation that 
is not a direct implementation of a physical quenching mech- 
anism. We leave a full comparison of the ezw model to a wide 
range of observables for future work, and focus here on the 
Hi properties of galaxies in our four wind simulations. 



3 THE HI MASS FUNCTION 

The most basic counting statistic for characterising the Hi 
content of galaxies is the Hi mass function (HIMF). Improv- 
ing 21cm observations have enabled the HIMF to be probed 
down to Hi masses appr oaching lO' Mq (e.g. IZwaan et ahl 
I2OO5I : iHavnes et allboill) . However, the redshift evolution 
of the HIMF remains poorly characterised, awaiting the 
next generation of facilities. Here we compare our simulated 
HIMFs with observations, to understand what constraints 
can be placed on wind models and how the HIMF is ex- 
pected to evolve. 
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Figure 3. Top: A comparison of the HIMF at z = in our four 
wind simulations — ezw (green), vzw (cyan), cw (red), and nw 
(blue) . The vertical dotted line is the stellar mass resolution limit 
of 1.4 X 10* Mq; this is roughly the H I mass resolution limit as well 
since M« and M^/ are comparable at the smallest masses (see 
Figure|2]l. The dashed li ne shows the HIMF from the o.40 sample 
of the ALFALFA survey llHavnes et al .I2OIII ') for comparison. The 
ezw and vzw model provide reasonable fits to the data. Bottom: 
Evolution of the HI mass function from z = 4 — > in the ezw 
model. The low mass end slope increases substantially at high 
redshifts. 



Figure |3l top panel, shows the HIMF for our four wind 
models. The vertical dotted line shows the stellar mass reso- 
lution limit of our fiducial r32n512ezw simulation. We have 
no formal Hi mass resolution limit, but since, as we will 
show in Figure 21 the Hi and stellar masses are similar for 
these small galaxies, this can roughly be considered as the 
H I mass resolution limit as well. Error bars (shown for ezw) 
depict cosmic variance as computed from the error of the 
mean HIMF in each of eight octants of the simulation vol- 
ume. 

For comparison, w e show the observatio nal determina- 
tion of the HIMF from iHavnes et al] (|201ll ) (the first 40% 
of the ALFALFA Survey; dashed line). We note that the 
ALFALFA observations select by HI mass, while we are ef- 
fectively selecting by stellar mass; however, the sensitivity of 
ALFALFA is well below what we can resolve, and hence un- 



less there is a dominant population of small Hl-free objects 
(which we will show later does not occur in our models), this 
comparison should be robust. 

Concentrating on the low mass end of the HIMF, we see 
that both the ezw and vzw models do a good job of matching 
the observed low mass end slope of a « — 1.34±0.02 (for th e 
"whole Q.40" sample, from Table 6 of iHavnes et al.ll201ll ). 
In detail, the predicted best-fit Schechter function slope for 
the ezw model is a = —1.31, while for the vzw model it 
is a = —1.45 and, therefore, the ezw model provides a 
slightly better fit. The massive end is relatively unaffected 
by quenching, as seen from the minor differences between 
the ezw model that includes quenching and the vzw that 
does not include it. 

Particularly remarkable is that the amplitudes of the 
HIMF in these two models are reasonably close to that 
observed. They are slightly above the observations at all 
masses, particularly at higher masses. We could in principle 
tune our values of A^Hi.iim and/or shrink the radial extent 
to which we associate Hi with a galaxy within reasonable 
uncertainties to match the data better, but we prefer to use 
well-motivated values for these quantities. We have in fact 
tried several other reasonable values for these quantities, and 
the net effect is generally to scale the HIMF in amplitude 
without changing the shape significantly. 

The agreement with the low mass end slope of the 
HIMF, while simultaneously matching the low mass end 
of the galactic stellar mass function, is a stringent con- 
straint that alm ost all galaxy formation mode l s hav e diffi- 
culty matching (IMo et al Il2"005l : ILu et al.ll2012l . |2013| ). Fun- 
damentally, in most models this arises because the low mass 
end slope of the dark matter mass function is quite steep, 
which is exacerbated by low-mass galaxies being more Hl- 
rich. In our simulations, the stronger outflows from low- 
mass galaxies strongly suppress the overall baryon content 
of galaxies. Even more critical is the impac t of w ind re- 
cycling. As described in lOppenheimer et al.l (|2010l ). wind 
recycling is preferentially stronger in high-mass galaxies be- 
cause the denser surrounding gas slows outflows more effec- 
tively via ram pressure. This results in a higher fraction of 
ejected gas being recycled back into higher-mass galaxies, 
there by yielding more sta r form ation in more massive sys- 
tems. lOppenheimer et al.l (|2010l ) showed that this flattens 
the low mass end slope of the stellar mass function (see Fig- 
ure [21 , and here we see that this effect is also important 
for the low mass end of the Hi mass function. We will dis- 
cuss further the differences between previous semi-analytic 
galaxy formation model results and our results in ^ 

Both the constant wind (cw) and no wind (nw) mod- 
els provide significantly poorer fits to the observations than 
ezw or vzw. Without winds, there is a surplus of Iow-Mhi 
galaxies, and a deficit of high-Mni systems. This is charac- 
teristic of th e HIMF in ma ny semi-analytic galaxy forma- 
tion models (|Lu et al.ll2012l ). The former discrepancy arises 
because there are no outflows to eject material from low- 
mass galaxies, and hence the b aryon fraction in these sys- 
tems is very large (jPavel [20091 ). in disagreement with ob- 
servations that show a small baryon content in dwarfs (e.g. 
iMcGaugh et al.ll2010l ). The latter arises because, without 
winds, gas becomes very dense in high-mass galaxies, which 
means that most of the gas is in molecular form, and 
much of it has been converted into stars. As discussed in 
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iDave. Oppenheimer. fc FinlatoJ (|201lh . the stellar content 
of galaxies in the no-wind model grossly exceeds that ob- 
served, and Figure [5] showed that it overproduces the num- 
ber density of galaxies at all the stellar masses probed. 

For the constant wind case, the low mass end slope is 
quite steep (a = —1.64). Moreover, there is a characteristic 
"bump" in the mass function at Mhi ~ 10^'' Mq. This bump 
is reminiscent of a similar bump in the stellar mass func - 
tion in this model l|Dave. Oppenheimer. fc Finlator|l201l[ ). 
which arises because wind recycling rapidly becomes impor- 
tant around this mass scale as the (constant velocity) winds 
are no longer able to escape the galaxy halo's potential well. 
Such a feature, which is generic to wind models that assume 
a constant outflow speed, is not observed in the HIMF (nor 
in the stellar mass function). 

iDuffv et all (|2012| ) employed an outflow model that is 
similar to this constant wind case. They find that the low 
mass end slope is fairly flat down to A/hi ~ 10^ M0, and 
was too shallow compared to data. In a comparable range, 
our HIMF also appears fairly flat - formally, the best-fit 
low mass end slope ignoring points below that mass is q = 
— 1.12, although a Schechter function is a poor descriptor. 
Our higher-resolution simulation probes further down the 
mass function, which enables us to see the steep low mass 
end. Overal l, our r esults for this wind model agree well with 
iDuffv et all (|2012D in the overlapping mass range, and both 
show that the constant wind model fails to match the low 
mass end of the HIMF. 

The bottom panel shows the redshift evolution of 
the HIMF out to z=4, focusing on our ezw simula- 
tion. The low mass end slope becomes progressively 
steeper with redshift, mi micking the behaviour seen in 
the s tellar mass function (|Dave. Oppenheimer. fc Finlatoi] 
I2OIII ). The low mass end slopes at 2 = 1,2,3,4,5 are 
-1.54, -1.79, -1.82, -1.99, -2.11, respectively (z = 3, 5 are 
not shown). This arises because wind r ecycling becomes in- 
creas ingly effective to lower redshifts l|Oppenheimer et al.l 
bold ), sinc e the r ecycling time is roughly constant at ~ 
1-3 Gyr (|Oppenheim cr & Dave 2008) for M. ~ 10^" Mq 
galaxies, which is a small fraction of the Hubble time today 
but comparable to the Hubble time in the early Universe. 
Hence the HIMF steepens rapidly out to 2: ~ 2, and then the 
steepening becomes more gradual, since at early times the 
effect of wind recycling, which is responsible for flattening 
the HIMF, is reduced. Meanwhile at the massive end, there 
are fewer galaxies at high redshifts simply because of the 
hierarchical nature of galaxy assembly. 

In summary, the HIMF provides a strong constraint on 
outflow models. The agreement of the ezw model predictions 
with the latest observed HIMF from ALFALFA represents 
a non-trivial success that has not previously been attained 
in hierarchical models of galaxy formation. The vzw model 
fares slightly worse but may still be within the overall un- 
certainties, while the constant wind and no wind cases fare 
poorly against data. More broadly, this indicates that in- 
cluding a well-motivated model for galactic outflows enables 
hierarchical structure formation models to produce an HIMF 
that is in very good agreement with data down to fairly low 
Hi masses. 
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Figure 4. Hi richness (= A/m/Af,) in galaxies from our four 
wind models: ezw (upper left), vzw (middle left), cw (upper right), 
and nw (middle right). Points show individual galaxies at z = 0, 
with central galaxies shown in blue and satellites shown in red. A 
running median for all galaxies is shown in green, while medians 
for the centrals and satellites are shown in blue and red, respec- 
tively. Circles show Hi fractions observe d in the GASS survey 
of nearby galaxies jCatinella et al.l I2OICII), squares sho w results 
from the Herschel Redshift Survey ijCortese et al.ll2011^ . and the 
dashed line shows an approxim ate fit to the locus traced by the 
Hl-selected ALFALFA survey l|Huang et al.ll2012f l. The bottom 
left panel shows a comparison of all four wind models with obser- 
vations at 2 = 0, while the bottom right panel shows the evolution 
of the median H I richness with redshift out to 2 = 2 in the ezw 
model. 



4 HI RICHNESS 

The H I richness, i.e. the H I mass relative to the stellar mass, 
provides a complementary characterisation of the H I content 
of galaxies. As we argued earlier, both the Hi and stellar 
masses are reasonably robust predictions of our models (un- 
like the molecular gas content), and hence the Hi richness 
should provide a meaningful discriminant between models. 

Figure |4] shows the Hi richness in our simulated galax- 
ies, as a function of stellar mass. The top four panels show 
our four wind models: ezw (upper left), vzw (middle left), 
constant (upper right), and no winds (middle right). The 
points in each panel show individual galaxies at 2: = 0, with 
blue depicting central galaxies and red depicting satellites. 
Galaxies with zero H I content are plotted along the bottom 



© 0000 RAS, MNRAS 000, 000-000 



HI in simulated galaxies 9 



of each panel; they are almost exclusively satellite galax- 
ies, which we will examine further in SjS] A binned median 
of log Mhi /M, for all galaxies is shown as the green line. 
We show errorbars on the median corresponding to the la 
spread for the galaxies within each mass bin. We also sep- 
arately show binned medians for the central and satellite 
populations. The vertical dotted line shows our galaxy stel- 
lar mass resolution limit. 

For comparison, mean Mhi/M* observations from the 
GALEX Are cibo SPSS Survey ( GASS) are shown as the 
large circles (CatincUa et al.lboij '). GASS has measured Hi 
in a stellar mass-limited sample down to M* ~ 10^°Mq. At 
lower masses, we show the results from the Herschel Refer- 
ence Survey (HRS) by ICortese et al] l|201ll . open squares), 
which uses literature Hi data but is approximately stellar 
mass-complete. We also show the fit to results from the AL- 
FALFA survey bv lHuang et al.l (|2012l . dashed line); since this 
is H I- mass selected, it is biased towards more H l-rich gala:x- 
ies, as is evident when compared to the M*-selected data. 
Our simulation results are most straightforwardly compara- 
ble to Mt-limited samples. 

All wind models broadly predict that Hi richness is 
anti-correlated with M*, as observed. However, in detail, the 
models show distinct differences; for clarity, just the overall 
medians are plotted in the lower left panel. Our ezw model 
produces a good agreement with the stellar mass-selected 
observations down to the lowest probed masses. The trend 
for the vzw model follows ezw, but shows slightly lowered 
H I richness particularly in smaller galaxies, indicating that 
more mass-loaded galactic outflows in low-mass galaxies are 
favoured. The no- wind case follows the trend of the data, but 
is too low by a factor of a few in Hi richness, showing that 
there is overly efficient conversion of gas into stars in this 
model. Finally, the constant wind model produces a distinct 
feature in Hi richness, mimicking the feature seen in the 
HIMF, owing to wind recycling. This agrees poorly with the 
observations, which has no comparable feature. Note that 
more highly mass-loaded outflows yield higher H I richness, 
not lower as one might naively expect, in part because such 
outflows suppress M, . 

The redshift evolution of Hi richness is shown in the 
lower right panel of Figure [4] for our ezw simulation. De- 
spite the rapid evolution in the HIMF, there is remark- 
ably little evolution in the Hi fraction at a given stellar 
mass from 2 = 0^2. The slow evolution in H I con- 
tent is a contrast to the much more rapid evolution ob- 
served in the molecular gas fraction (e.g. lTacconi et al]|2010l : 
iGeach et al. I I2OII ), although some of that evolution may 
reflect uncertainties in assessing the molecular gas content 
from CO emission (e.g. Narayanan. Bothwell. fc Davell20l3 : 



iBolatto. Wolflre fc Lerovll2013l '). Also, the lack of evolution 
in Hi richness indicates that the evolution in the HIMF 
discussed in i|3] mostly reflects the evolution in the stellar 
and/or halo mass functions. We will see in fj7]that the lack 
of evolution in global H I content is also a generic prediction 
of our models. 



5 H I IN SATELLITE GALAXIES 

The Hi content of galaxies is seen to vary substantially 
with environment. In the densest environs such as clus- 
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Figure 5. The top panel shows Hi richness (= Mhi/M*) for 
satellite galaxies versus halo mass at 2; = 0, for our ezw model. 
The satellites are split into low (blue), intermediate (green), and 
high (red) stellar mass bins, as indicated. Satellites with no H I 
are plotted at —3, with an artificial spread to aid visualisation. 
Lines show a running median of only the galaxies that have H I 
(i.e. ignoring the points plotted along the bottom). The bottom 
panel plots the fraction of Hl-free satellites within those same 
mass bins, colour coded as above. Solid lines show results from 
the ezw model that includes an ad hoc quenching model at high 
halo masses, while dotted lines show the equivalent plots for the 
vzw model that does not include such quenching. 



ters, ram pressure can remove Hi from infalling galaxies, 
and it i s long known that cluster galaxies are deflcient in 
Hi (e.g. iHavnes et al.lll984 ). Even in less extreme environ- 
ments where ram pressure stripping is expected to play less 
of a role, the Hi content appears to be anti-corr elated with 
local density (|Robertson. Shields, fc Blandl2012l ') ; other pro- 
cesses such as strangulation or harassment may also be play- 
ing a signiflcant role. In general, such environmental pro- 
cesses are expected to preferentially lower the Hi content 
of satellite galaxies, particularly those within larger halos 
that are expected to have hotter ambient gas. In this sec- 
tion we examine how the H I content of satellites varies with 
halo mass, which can be considered as a rough proxy for 
environment. 

Returning to Figure [l] it is clear that satellites (red 
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points and curve) have a lower Hi richness than central 
galaxies (blue) of the same stellar mass. At the high mass 
end, the typical difference in median Hi richness between 
central and satellites is a factor of ~ 2 — 3, in all the wind 
models. Below M* < IQ'^-^Mq, the satellite Hi fraction 
drops very quickly, as there are numerous low-mass H I- free 
satellites in our simulations (plotted along the bottoms of 
the panels). While the majority of galajdes of low ma ss are 
still centrals (|Dave. Qppenheimer. fc Finlatoij l201ll ). and 
hence the overall median tracks that of central galaxies, it 
is clear that low-mass satellites in particular are highly de- 
ficient in Hi for their mass. 

Figure [S] examines these trends more closely. The top 
panel shows the H I richness in satellite galaxies in our ezw 
simulation at z = 0, divided into three stellar mass bins of 
1.4 X 10* - 10^ M©, 10^ - 10^° Mq, and > IQ^°Mq. Galaxies 
with H I richness less than 10"'^ (virtually all of which are 
Hl-free) are plotted along the bottom at —3; lines show a 
running median not including these H l-free galaxies. 

This plot shows several key trends. First, at all masses, 
the satellites that have H I show a mild trend of being more 
Hl-rich in lower mass ha los. Since the halo mass traces 
stella r mass in our models (|Dave. Qppenheimer. k, Finlatoi] 
this basically reflects the fact that lower stellar mass 
galaxies have higher Hi richness (as seen in Figure [4]). Sim- 
ilarly, lower- mass satellites are more Hi rich, again follow- 
ing the trend for centrals discussed in the previous section. 
Overall these trends appear to reflect the H I content of the 
satellite galaxies when they were still centrals. 

However, a clear difference between central and satellite 
galaxies is that there are many more satellites that are de- 
void of Hi, particularly at higher halo masses, i.e. the points 
along the bottom of the plot appear much more frequently in 
the more massive halos. In fact, the distribution of H I rich- 
ness in satellites appears to be bimodal. Even in large halos, 
there are still some satellites that have substantial Hi, and 
their Hi content is not grossly different than that of satel- 
lites in smaller halos. However, the fraction of satellites that 
are Hi poor increases sharply in higher mass halos. 

To quantify this, we plot in the bottom panel of Fig- 
ure [5] the fraction of galaxies with Hi richness < 10^'^, as 
a function of halo mass, divided again into the three satel- 
lite galaxy stellar mass bins as in the top panel. Here, we 
clearly see that at any halo mass, lower-mass satellites are 
more likely to have had their Hi content strongly reduced. 
Moreover, this trend is a very strong function of halo mass, 
with satellites of all masses being much more likely to be 
H l-poor if they lie within a more massive halo. 

The strong bimodality suggests that the process that 
renders satellites Hi poor happens on a relatively short 
timescale compared to the infall timescale into the halo 
(which is roughly comparable to the halo's dynamical time 
of several Gyr). In future work we plan to investigate the 
detailed dynamical processes that remove the H I from satel- 
lites in our simulations. 

The dotted lines in the bottom panel of Figure [5] show 
the analogous results for our vzw simulation. In massive ha- 
los where the fraction of H i-poor satellites is substantial, the 
outflow model is identical between vzw and ezw simulations, 
but the ezw simulation includes our quenching prescription. 
The fact that the dotted and solid lines are very similar 
indicate that the trends seen in the Hi content of satel- 



lites are not being set by our ad hoc quenching prescription, 
but more likely by the fact that halo masses above 10^^ Mq 
tend to contain much more hot halo gas th at can more 
strongly impact satellite s moving through it l|Keres et al.l 
I2OO5I : iGabor et all I2OIII ) . Recall that our quenching pre- 
scription applies only to galaxies with a high velocity dis- 
persion, which is typically only the central galaxy in the 
halo, and hence the only direct impact on smaller satellites 
would be from the extra heat being added to the halo gas; 
evidently this has a minimal effect on the satellites in our 
simulations. 

What about the central galaxies? In general, very few 
of the central galaxies are devoid of H I, as seen in Figure [l] 
However, although we don't highlight it, there exists a small 
population of lower- mass (A/« < 10^" M©) cent rals that are 
H i-po or. This is related to what was seen bv lGabor et al.l 
\2Q1% . who found numerous low-mass central galaxies on 
the red sequence. These turned out to be galaxies that re- 
side just outside, i.e. within several virial radii, of more mas- 
sive halos. In the spherical overdensity algorithm we use to 
identify halos, such galaxies are identified as centrals, al- 
though the influence of the larger gala xy's halo can exten d 
to well beyond its virial radius (e.g. iMoldar et al.l [2OO9I ). 
Hence these H l-poor centrals could be impacted by the ex- 
tended environment of a nearby larger halo, or else they 
could be former satellites whose orbit has taken them out- 
side the nominal virial radius. 

In summary, halo mass plays an increasingly important 
role in setting the Hi content of satellite galaxies. This is 
particularly seen by the strongly increasing fraction of H l- 
poor satellites as a function of halo mass. At a given halo 
mass, low-mass satellites have a greater chance of having 
their Hi removed. Satellites that have not had most of their 
H I removed lie along similar relations to satellites in lower- 
mass halos. These results suggest that the process by which 
Hi is removed from satellites in our simulations acts fairly 
quickly, and preferentially on smaller galaxies. Comparing 
these predictions to observations can help constrain such 
Hi removal mechanisms. 



6 H I DEFICIENCY 

The Hi richness of galaxies is observed to correlate with 
a variety of galaxy properties (besides stellar mass). In 
the previous section we showed that satellite galaxies in 
higher mass halos were increasingly stripped of their Hi. 
But even galaxies that still have substantial Hi show 
correlations of their Hi content with properties such 
as environme nt (e.g. 



Cortese et al. ^ 201 ll ). and metallic- 
ity (Robertso n. Shields, fc Blanc 2012,). Here we examine 
the second-parameter trends of Hi with environment, star 
formation rate, and metallicity in our simulations to pro- 
vide insights into the physical drivers that establish the H I 
content of galaxies. 

To set a theoretical context for this discussion, we re- 
ca ll the equilibrium model for galaxy gr owth as described 
in lPave. Finlator. fc Oppenheimen l|2012l 'l. which presents a 
simple physical scenario that can account for the relation- 
ships between key galaxy properties, as well as the scatter 
around those relationships. In the equilibrium model, accre- 
tion onto a galaxy from the cosmic web is fairly quickly 
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Figure 6. Left panels: Specific SFR, metallicity, and environment 
vs. stellar mass (top to bottom) from our ezw simulation at ^ = 0, 
with the blue and red points indicating galaxies that have above 
and below (respectively) the median Hi richness at that given 
Mf. The green line shows the running median. The cyan line 
shows the median for the vzw model, although the points are 
not shown. The black solid line in the middle-left panel plots 
the median observed relation from l lTremonti et al.ll2004) with 
the dashed lines indicating the 95% contours. Hl-rich galaxies 
at a given M« tend to have higher sSFR, lower metallicity, and 
live in less dense regions. Right panels: The deviation from the 
median sSFR, metallicity, and environment (top to bottom) at 
a given vs. the deviation from the median H I richness. The 
magenta lines show power law fits to the deviations, having slopes 
of 0.31, —0.26, and — 0.56, respectively. 



processed into either stars or an outflow, resulting in a 
slowly-evolving gas reservoir. This results in fairly tight re- 
lations between stellar mass and star formation rate (the 
so-call ed main sequence , e.g. IPavel 120081 ) . gas-phase metal- 
licity jFinlator_&_^W3j200^)j_^^ gas con- 
tent (|Dave. Finlator. fc Oppenheimeij|201ll) . The predicted 
relations are in good agreement with the hydrodynamic sim- 
ulation results. 

Stochasticity in the inflow rate causes scatter 
about these equilibrium relations. As described in 
iDave. Finlator. fc Qppenheimerl (|201ll ). an accretion event 
such as a minor merger will cause an increase in gas con- 
tent that raises the star formation rate, while simulta- 
neously lowering the metallicity. Conversely, a lull in ac- 



cretion, or diminished accretion owing to a galaxy be- 
coming a satellite in a larger halo, will cause the ex- 
isting gas to be consumed, resulting in a lower gas 
content, a lower SFR, and a higher metallicity. Hence, 
for example, departures from the mass-metallicity rela- 
tion are expected to be inversel y correlated with SFR . 
Such a trend has be en observed l|Lara-L6pez et al] |2010| : 
iMannucci et al.l |2010| ). and is known as the fundamen- 
tal met allicity relation; our simulations na turally predict 
this (Da ve. Finlator. fc Oppenheimej l201lh . Our models 
analogously predict that deviations in the star-forming (i.e. 
mostly molecular) gas con tent will inversely correlate with 
devia tions in metallicity (jPave. Finlator. fc Qppenheimej 

(ml). 

We now extend this to consider the impact of inflow 
stochasticity on the H I reservoir of galaxies. The left panels 
of Figure |6] show the specific SFR, the metallicity, and the 
local galaxy density averaged over 1 Mpc spheres as a func- 
tion of stellar mass in our ezw simulation at 2 = 0. Metallic- 
ity here is computed as the SFR-weighted oxygen metallic- 
ity |Finlator fc Davil2008l : iDave. Finlator. fc Qppenheimej 
boiC " The observed mass-metallicity relation (MZR) from 
SDSS (|Tremonti et a l. 2004) is shown as the solid black line 
with enclosing 95% contours indicated by the dashed lines. 
We have converted these metallicitie s to solar units assum- 
ing a solar oxygen abundance from lAsplund et al.l (|2009l ) , 
namely 12 -I- log [0/H]g = 8.70. The green line shows a run- 
ning median in each panel. We also show as the magenta 
line the median for the vzw simulation (without showing 
the individual points), for comparison. 

The median sSFR drops slowly with Af* at at 
both small and large stellar masses, with a peak at 
M. - 5 X 10^ M0; there are minimal differences be- 
t ween the ezw and vzw models. The vzw line follows that 
in lDave. Oppenheimer. fc Finlatoil (|201ll ). who used lower- 
resolution simulations that only probed down to M, ^ 
10® Mq, while here the turnover is much more apparent. 
We note that SDSS obser vations do not indicate such a 
turnover (|Salim et al.ll2007l ). which reflects a generic prob- 
lem in hierarchical models that star for mation in dwarf s 
peaks too early and is too low today tWeinmann et aLll2012l ). 

The metallicity is tightly correlated and increases with 
M. untfl M. ~ 5 X 10^° M0. The ezw model produces 
a steeper mass-metallicity relatio n (MZR) than t he vzw 



mod el, for reasons discussed in Finlator fc Davel (|2008l ) 



and iDave. Finlator. fc OppenheimerT i 20121 ). The metallic- 



ity Z oc 1 /{I + rj) (in the absence of wind recycling), with 

7) oc l/(j oc Mt^^^ for all galaxies in the vzw model and 

l/cr^ oc M~^^^ for the ezw model at small masses. For 

1/3 

77 ^ 1 as at small masses, this toy model yields Z oc 
for the vzw model and Z oc Af^'^'^ for the ezw model. This 
represents an asymptotic slope, and since the ezw model fol- 
lows the vzw model scalings at high masses, the actual MZR 
in the ezw model is not as steep as this, but it is clearly 
steeper than in the vzw model. The vzw model appears to 
provide a better fit to the SDSS data, although metallicity 
measures are sufficie ntly uncertain that the ezw prediction 
cannot be ruled out (|Ellison et al.l [20081 ). Also, a recent de- 
termination of the M ZR fr om direct metallicity measures by 

Andrews & Martinil l|2013l ) suggests a steeper MZR slope at 

1/2 

small masses, closer to Z oc A/, . 
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The simulated galaxies in Figure [6] are colour-coded 
by H I richness: the blue points show galaxies that have a 
higher-than-median Hi richness at that stellar mass, and 
the red points show the converse. It is clear that Hi rich- 
ness has a strong second parameter correlation with all these 
properties: H l-rich galaxies at a given mass also tend to have 
a higher sSFR, a lower metallicity, and a lower local galaxy 
density. 

To reiterate the physical scenario: When accretion hap- 
pens, it results in galaxies with higher Hi richness than typ- 
ical. This high H I richness also corresponds to a lower-than- 
normal gas-phase metallicity. This suggests that the H I con- 
tent of galaxies can be an i ndicator of recent a ccretion, as 
argued from observations bv lMoran et al.l l|2012l ). who show 
that the low metallicities in Hl-rich systems most strongly 
appear at the outskirts of disks. 

We can quantify these trends using deviation plots, 
which we show in the corresponding right panels of Figure[6l 
This shows the difference between a galaxy's sSFR, metal- 
licity, and density with the median value of these quantities 
at a given M*, i.e. approximately the green line in the left 
panels (although in detail we employ a spline fit to produce 
more smoothly varying deviations), versus the deviation in 
Hi richness with its median value at that M, . Points in 
these panels are colour-coded by centrals (blue) and satel- 
lites (red). The green lines show power-law fits to these de- 
viations for the ezw model while the magenta lines show the 
same for the vzw model (though the individual points are 
not shown). 

These right panels quantify the second-parameter 
trends from the left panels: galaxies with high Hi content 
also have a high sSFR, a low metallicity, and are in low 
density environments. A power law provides a reasonable fit 
to the sSFR and metallicity data, with slopes of 0.3f and 
—0.26, respectively. The vzw model yields virtually identical 
slopes, even though the trends (such as the MZR) are notice- 
ably different. This occurs because the trend in the scatter 
arises from inflow stochasticity, which has little to do with 
the outflows. The central and satellite galaxies do not show 
significantly different trends in sSFR or metallicity. 

The metall icity dependence can be compared to re- 
cent work by [Robertson. Shields, fc Bland (|2012l ). who 
looked at the Hi deficiency parameter, DEF, relative 
to the deviation from the expected oxygen abundance 
for cluster and field galaxies, finding slopes of —0.25 ± 
0.f2 and —0.41 ± 0.14, respectively. Our simulated galax- 
ies are more comparable to field galaxies, which show 
a higher slope but are st ill with in their uncertainties. 
[Robertson. Shields, fc Bland (|201^') compared to results 



ronment deviation plot (bottom right) markedly separates 
the central and satellite galaxies. This environmental de- 
pendence arises because the environment can impact the 
inflow rate into satellites. High density regions associated 
with halos with masses > 10^^ Mr will conta i n substantial 

2009; Gabor et all 



from the simulations in Dave, Finlator. fc Oppenheimej 
l|201ll ). which yielded a similar slope. Although the defi- 
nition of DEF differs from the deviation we plot here, the 
qualitative agreement in the trend indicates that the basic 
physical model of a slowly-evolving equilibrium in gas con- 
tent and metallicity appears to be broadly consistent with 
the observations. 

Finally, we consider the environment in the bottom two 
panels of Figure [G] The median local galaxy density is in- 
dependent of mass until the very highest mass systems in 
our volume at M* > 5 x 10^° M©. Unlike in the sSFR and 
metallicity deviation plots, where the centrals and satel- 
lites do not significantly deviate from one other, the envi- 



amounts of hot gas (jKeres et al.r 2005. _ . 

l 2012i') that can retard accretion (jDekel fc Birnboim 20061 : 
iDave. Finlator. fc Qppenheimerl |2012|) . Hence galaxies in 
such regions, particularly satellites, will have lower accre- 
tion rates compared to field galaxies at the same mass, and 
hence less Hi. Overall, the full galaxy sample trends to- 
wards having less gas-rich galaxies in denser environments, 
but this is driven almost entirely by the satellites. Fitting 
a power-law formally yields a slope of —0.56 (shown as the 
magenta line), though a power law does not appear to be a 
particularly good fit. 

In summary, the correlated deviations between sSFR, 
metallicity, and gas content provide quantitative constraints 
in the way in which galaxies oscillate about their equilibrium 
relations owing to stochastic accretion. The resulting devia- 
tion slope reflects the correlation between the instantaneous 
gas inflow rate, its conversion to neutral hydrogen, and the 
infalling metallicity, as mediated by the local environment. 
Such constraints provide valuable discriminants between dif- 
ferent physical scenarios for galaxy growth. The consistency 
with available observations, albeit preliminary and qualita- 
tive, indicates that the stochasticity in the inflow as expected 
from cosmological accretion is able to broadly explain the 
second parameter trends observed in the relationship be- 
tween Hi content and other physical galaxy properties. 



7 Qm EVOLUTION 

Most neutral hydrogen in the cosmos exists in and around 
galaxies, between the regime where self-shielding happens 
in the outskirts of disks, to where gas becomes molecular- 
dominated within star-forming regions. While it is currently 
infeasible to probe the evolution of Hi in galaxies via 21cm 
emission to high redshift^, other avenues have been em- 
ployed to measure cosmic Hi evolution, such as the abun- 
dance of damped Lya systems (DLAs). In this section we 
examine the evolution of the cosmic H I density in our sim- 
ulations directly from the galaxy population, and compare 
it with various observational determinations. 

Figure [7] shows the evolution of JIhi from z = 5 — )■ 
in our vzw model, for all galaxies above our approximate 
Hi mass resolution limit of Mm > 1.4 x lO^M© (solid 
blue line), as well as for galajcies with Mhi > W^Mq 
(dashed) and Mhi > IO'^'^Mq (dotted). For comparison, 
we also show the evolution of the stellar mass density SI, 
as th e red line. Observat i ons o f flm &t z > 2 are shown 
from iProchaska fc Wolfd (|2009l ) derived from DLA abun- 
dances (circles with errorbars). At 0.5 < z < 2, we show a 
box that encompasses the Icr range from ^ Rao et al.l (|2006l ). 
who used strong Mgll systems as a proxy for DLAs. Fi- 
nally, th e triangle at z = shows SIhi measured from AL- 
FALFA (jHavnes et al.lfioill '). Note that our prescription for 



^ This is a key g oal of the upcoming L ADUMA survey using the 
McerKAT array l lHolwerda et al.ll201lh 



© 0000 RAS, MNRAS 000, 000-000 



HI in simulated galaxies 13 



-3 

a -4 
-5 
-6 




o 



n 
n 
n 



HI,>1.4e8 
HI,>le9 



HI,>lelC 



Hl.tol 

n. 



1 2 3 

Redshift 



4 



Figure 7. QhIi the fraction of cosmic mass density in Hi, from 
2 = 5 — > from our ezw simulation. The blue solid line shows 
Ohi from all resolved galaxies (M, > 1.4 X 10* Mq) the blue 
dashed and dotted lines show Qhi in galaxies with M« > 10® Mq 
and Mt > 10^^ Mq, respectively. The cyan line shows Ohi from 
extrapolating a Schechter function fit to the HIMF at each red- 
shift down to Mhi = IO^Mq. Magenta line shows the total 
H I mass density in the entire volume. For comparison, the red 
line shows the mass density in stars in this model. Observa- 
tions are shown at 2 < z < 4.5 (circles with errorbars) from 
iProchaska fc W olfe (2009|), at 0. 5 < 2 < 2 (box encompassmg 
la limits) from Rao et all ||2006|) . and at 2 = (triangle) from 
iHavnes et al] ll201lfK 



determining self-shielding calibrated to radiative line trans- 
fer simulations as described in ^does reasonably well pre- 
dicting nei at 2: = 0; we do n ot adjust any parame ters to 
match this data as was done in lPopping et all l|2009t ). 

Our simulations predict that Qhi is remarkably con- 
stant from z = 5 — >■ 0, essentially unchanging to within 50% 
over these 12 billion years. This is in contrast to the dra- 
matic increase of SI* over that time interval; over 98% of 
all stars in this model form since 2 = 4, with more than 
80% since z = 2 (in b road agreeme nt with recent obser- 
vational estimates; e.g. iLeitnerl |2012| ). It is also quite dif- 
ferent than the r apid evolution in cosmic s tar formation 



rate density (e.g. iHopkins fc BeacomI 20061: Fardal et al. I 



200'itl and molecular gas fraction " iTacconi et al.l l2010l : 



Geach et al1l201ll : iTacconi et al]|2012l ). which both increase 
by an order of magnitude or so out to z = 2 (though 
see [Narayanan. Bothwell. fc Pavel |2012| . who argue for less 
growth in the observations owing to variations in the CO- 
to-H2 conversion factor). This emphasises that Hi repre- 
sents a transient reservoir in the journey of gas from the 
ionised IGM to stars forming deep within galaxies, and 
is not directly proportional to the amount of stars being 
formed or molecular gas present. This cautions against over- 
interpreting quantities such as the "Hi star formation effi- 
ciency", i.e. the SFR divided by the Hi mass, since the two 
quantities are not directly related. 

Subdividing SIhi into different Hi mass bins, we see 
that at z = 0, most of the cosmic Hi is in galaxies with 
10'^ < Mhi < 1O^°M0. Going back in time, Hi shifts to- 
wards lower mass galaxies; hj z — 2, only half the H I 
is in Mhi > lO'A/0 systems. In contrast, only a very 



small portion of the cosmic H I ever resides in galaxies with 
Mm > 10^" Mq, since higher mass galaxies quickly drop off 
in their H I richness (Figure^). Thus quenching mechanisms, 
which primarily affect high-mass galaxies (e.g. iGabor et al.l 
l2012h . are expected to have almost no effect on Qm . We also 
show, as the magenta line, the total Hi mass density in the 
entire simulation volume. The difference between the solid 
magenta and blue lines reflects H l that is not within resolved 
galaxies (including the diffuse IGM). This extra contribu- 
tion is a couple percent at low redshift but rises to ~ 50% 
at 2 = 5, which again reflects the increasing contribution of 
small (unresolved) galaxies to the global Hi budget at high 
redshifts. 

Observations of Qm likewise indicate very little evolu- 
tion from 2 ~ 4 ^ 0. Different observational tracers gen- 
erally agree to within a factor of approximately 2, and to- 
gether indicate essentially no change in Qui for the past 
12 billion years, in broad agreement with our predictions. 
Our simulation begins to underpredict Qui at z > 3 when 
compared with the DLA data. The discrepancy could be 
physical or numerical. A physical explanation could be that 
lower metallicity galaxies at early epochs will tend to have 
less efficient conversion of their atomic gas to molecular gas 
owing to lower cooling rates and harder interstellar radia- 
tion, and hence our loc ally-calibrated prescription for i?moi 
from iLerov et al.l l|2008l ) may not be appropriate. 

A numerical explanation for the discrepancy may be 
that we do not fully resolve galaxies with Af, < 1.4x 10* M©, 
and there may be substantial contributions to Qui from the 
very smallest galaxies. At low-z, the low mass end of the 
HIMF is fairly shallow, so the expected contribution from 
lower mass galaxies is small. But at higher redshifts, Figure[3] 
shows that the slope becomes substantially steeper, meaning 
that the additional contribution from unresolved systems 
could be large. 

We can crudely correct for this by fitting a Schechter 
function to the HIMF at each redshift as we did in SjS] and 
then integrating down to some chosen lower mass limit. As 
an illustration, we integrate down to A4ui = IO^Mq, which 
is around the lowest mass observable at \ow-z in large sur- 
veys such as ALFALFA. The result of this Schechter fit ex- 
trapolation to Mhi = IO^Mq is shown as the cyan line in 
Figure [T] This results in a negligible correction at low-z, 
but at higher redshifts the correction can be up to a factor 
of three, which agrees better with the high-z DLA results. 
We caution that this exercise is intended to be illustrative, 
since our A/hi limit was chosen rather arbitrarily, and it may 
not correspond to the effective Hi masses probed by DLA 
systems. Furthermore, our Schechter function fits become 
increasingly uncertain at higher redshifts, owing to the lack 
of dynamic range in our simulations. Nonetheless, this illus- 
trates that plausible corrections down to lower H I masses 
can bring our predicted Qui into better agreement with the 
observations by preferentially increasing the high-redshift 
Hi mass density. Given the crudeness of our prescriptions 
for determining H I content, agreement at this level is en- 
couraging. 

In summary, our simulations generally predict a very 
slowly-evolving cosmic H I mass density, in broad agreement 
with the observations. However, we note that the HIMF 
actually evolves rather considerably (Figure [S]), with many 
more 1ow-A/hi galaxies and fewer high-Mni ones at high 2. 
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Figure 8. Comparison of tiie Hi mass function (top panel) and 
Hi richness (bottom panel) in our fiducial 32/i~^Mpc, 2 X 512^ 
particle simulation (green) compared to a 48/i~^Mpc, 2 X 384'^ 
particle simulation (red) with 8x poorer mass resolution. Res- 
olution convergence is good for H I richness, while the HIMF is 
increased in the higher-resolution simulation, but is still within 
the formal uncertainties. We show, as the magenta line in the top 
panel, the r48n384 HIMF where the Hi masses (not ■!>) have been 
multiplied by 1.5; this results in a nearly identical mass function 
to the r32n512 case. 



It is something of a coincidence that these two variations 
roughly cancel when summed to give the global Hi mass 
density. Nonetheless, this emphasises the transient nature 
of the H I reservoir around galaxies, which does not build up 
hierarchically in the way that the stellar mass does, but in- 
stead responds to evolution in the ionising background and 
the ISM pressure. 



8 RESOLUTION CONVERGENCE 

In this study we have employed some of the highest- 
resolution hydrodynamic simulations of random cosmolog- 
ical volumes ever evolved down to z = Q. Nonetheless, 
it is important to assess whether our resolution is suffi- 
cient to robustly predict the H l properties of galaxies. Here 
we examine our two most basic Hi statistics, namely the 



Hi mass function and the Hi richness vs. stellar mass, in 
a simulation that has 2 x 384^ particles in a 48/i'^Mpc 
volume (r48n384) , described in our previous work (e.g. 
Oppenheimer et al. 2010l: Dave. Qppenheimer. fc Finlatorl 



201ll : Dave. Finlator. fc Oppenheimerl 2011 ). This simula- 
tion has 8 times poorer mass resolution and 2 times poorer 
spatial resolution, albeit in a volume that is 3.4 times larger. 
We use the identical momentum-driven wind (vzw) and cool- 
ing model in both simulation^. 

Figure [8] shows the Hi mass function (bottom) and the 
Mhi/M, ratio vs. Af, (top) and in these two simulations, 
r32n512vzw in green and r48n384vzw in red. These are anal- 
ogous to Figures |4] and [S] respectively, and the observations 
are plotted as in those figures. 

For the Hi fraction, the resolution convergence is excel- 
lent; the running median lines for the high and low resolu- 
tion simulations lie on top of each other in the mass ranges 
that overlap. It shows that the anti-correlation between H I 
richness and stellar mass is converged even at significantly 
poorer mass resolution. Hence, at least over this modest dy- 
namic range in mass resolution, the H I richness of galaxies 
is a well-converged property in these simulations. 

In contrast, the HIMF (top panel) is less well converged. 
In general, the number density at a given H I mass increases 
with resolution. As shown by the magenta line, the two mass 
functions lie on top of each other if the Hi mass in each 
low-resolution galaxy is multiplied by a factor of 1.5. Yet 
the good agreement with the H I richness plot would suggest 
that the differences in the HIMF between these models ac- 
tually reflects differences in the stellar mass function. Indeed 
this is seen to be true (not shown), which may arise owing 
to increased wind recycling in higher-resolution simulations, 
or simply owing to the cosmic variance within the smaller 
volume; we leave this investigation for the future. Unfor- 
tunately, we do not have even higher resolution simulations 
within a sufficiently large volume to assess whether this con- 
verges with increased resolution. In any case, the differences 
in "I> between the models at any given Af, is smaller than 
the variance within that bin. 

Overall, the resolution convergence for H I properties is 
good but not ideal. For a given A/* the Hi content is very 
robust, but the number density of galaxies at a given Hi 
mass is less well converged. For that reason, ^hi is also 
slightly uncertain, as there is a 30% difference between 
the z = value predicted by the high and low resolution 
simulations. We do not yet know whether these differences 
reflect a resolution effect or a volume effect; we would need a 
wider suite of simulations to assess these issues. While these 
discrepancies are not trivial, they do not signiflcantly affect 
the main conclusions of this paper. 



9 SUMMARY AND DISCUSSION 

The H I content of galaxies offers a unique glimpse into 
baryon cycling processes of inflow and outflow that are be- 
lieved to govern galaxy evolution. Hydrodynamic simula- 
tions can provide a way to connect H I to such baryon cycling 



^ We use the vzw wind model because we already have the 
r48n384vzw simulation in hand, and there are usually only minor 
differences between vzw and our now-favoured ezw model. 
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processes, as well as to help interpret observations within a 
hierarchical galaxy formation context. In this work, we have 
examined the Hi properties of galaxies in cosmological hy- 
drodynamic simulations with four different outflow prescrip- 
tions. Our key conclusions can be summarised as follows: 

• The Hi mass function is strongly affected by the in- 
clusion of galactic outflows. Without outflows, our simu- 
lations yield too many galaxies at small Hi masses, and 
too few at large Hi masses. Our model with a constant 
wind speed (cw) introduces a feature in the HIMF around 
Mhi ~ Mq that is not seen in the observations. Our 
models with momentum-driven wind scalings produce an 
HIMF that broadly agrees with the observations, partic- 
ularly at the low mass end down to our mass resolution 
limit of 10^ Mq. Switching to energy-driven scalings at 
low masses (crgai < 75 km/s) further improves the agree- 
ment, and generates a low mass end HIMF slope of —1.3. 
The HIMF is slightly above the data at all masses, but likely 
within the systematic uncertainties regarding how we com- 
pute Hi masses for our galaxies. The low mass end slope 
becomes progressively steeper with redshift out to z ~ 3, 
then remains constant at around —2.1 at higher redshifts. 
As an aside, we showed that the ezw model also provides 
an excellent match, within the quoted uncertainties, to the 
galaxy stellar mass function. 

• All models show that low-mass galaxies are more H I 
rich, relative to their stellar content, than high mass galax- 
ies. The exact shape and amplitude of the relation de- 
pends slightly on the wind model. All the models broadly 
agree with the GASS data that probes down to M, ~ 
10^° Mq, but the different wind models predict different 
trends to lower masses. The constant wind model predicts 
significantly lower gas richness than our vzw (momentum- 
driven wind scalings) or ezw models at M, < 10^'^ M©, re- 
flecting a_strong_su£2I6Ssion^ inflow owing to energetic 
winds (Ivan de Voort et al. '2010!}_ that lik ewise suppresses 
specific SFRs l|Dave, Op pcnhci mer. fc Finl ator 2011). Cur- 
rent observations favour a continuing increase of gas richness 
to lower masses, and our ezw model again fares slightly bet- 
ter than vzw at low masses, and both match much better 
than the cw model. There is essentially no redshift evolution 
in the H I richness as a function of stellar mass, showing that 
the evolution of the HIMF is interconnected with the evolu- 
tion of the galaxy stellar mass function. 

• Galaxies with a high Hi content tend to have lower 
gas-phase metallicities and higher star formation rates at 
a given M*. The deviation in sSFR and metallicity ver- 
sus the deviation in H I richness can be fit by power laws 
with slopes of 0.31 and —0.26, respectively. The latter slope 
compares favourably with existing observations of the H I 
deficiency parameter versus metallicity deviation. These de- 
viation trends can be understood within the context of an 
equilibrium model in which galaxies tend to have equilib- 
rium relations in star formation rate, metallicity, and gas 
content, perturbed by stochastic accretion events that result 
in correlated deviations in SFR, metallicity, star-forming gas 
mass, and Hi mass. This suggests that Hi content is an in- 
dicator of recent accretion events in the outskirts of galaxies 
that stimulates star formation. 

• Environment plays an important role in governing H I 
content, causing an increasing suppression of the Hi content 



of satellites in more massive halos. The suppression is rapid, 
as indicated by the bimodal distribution of the H I fraction 
at a given halo mass. Lower mass satellites are more likely 
to have their Hi removed. The deviation of Hi with envi- 
ronment is anti-correlated with the deviation in Hi richness, 
driven by H l-poor satellite galaxies in denser regions. 

• The global H I mass density evolves slowly from 2: ~ 5 — >■ 
0, in broad agreement with the observations of DLAs and 
other measures. This is in contrast to the stellar mass den- 
sity, which grows substantially over that interval, highlight- 
ing the transient nature of the H I reservoir around galaxies. 
Today, the majority of the cosmic Hi mass is in galaxies 
with 10^ < Mhi < IQ^AfQ. 

• We briefiy examine numerical resolution convergence 
over a modest span of a factor of 8 in mass resolution be- 
tween two simulations with momentum-driven winds. The 
H I richness shows better convergence than the H I mass func- 
tion, which increases as the resolution is improved, equiva- 
lent to an increase in H I mass by a factor of 1.5. Hence reso- 
lution convergence has not been demonstrably achieved yet, 
but the differences remain comparable to or smaller than 
the statistical uncertainties. 

Simultaneously matching the observed low mass end 
of the galactic stellar mass function and the low mass end 
slope of the HIMF has remained elusive in simulations and 
even in semi-analytic models (SAMs) of galaxy formation 
where o ne has much m o re freedom to adjust the physical 
model (|Mo et alllioosi : ILu et"alll2012l . I2OI3I ). Fundamen- 
tally, in most models this arises because the low mass end 
slope of the dark matter mass function is quite steep, which 
is exacerbated by low- mass galaxies being more Hl-rich. 
iMo et al I l|2005h argue that it is not possible to match the 
low mass end of the HIMF in a very broad range of CDM 
based galaxy formation models if one only makes two ba- 
sic assumptions: First, that the original gas distribution has 
the same specific angular momentum as the dark matter 
ha lo, which it conserves as i t forms an exponential disk (as 
m [MoZ Mao, fc White! 1 19981 ): second, that stars only form 
in gas above a c ritical surface density Y^crit, as observed 
(|Kennicuttl Il989l ) . Almost all SAMs make these same as- 
sumptions. They populate some fraction of the cold dark 
matter halos with exponential gas disks and assume that 
any region originally above Ecrit has its surface density low- 
ered to Tiarit- Making the (conservative) assumption that 
half of the gas is neutral, they overproduce the observed 
HIMF at the lo w mass end by a factor of more than five. 
ILu et all 1I2OI2I ) confirm this result by adjusting all their 
SAM parameters using a Bayesian inference technique to 
match the observed low redshift _K'-band galaxy luminosity 
function, and find that they overpredict the observed HIMF 
by a large factor. 

ILu et al] (|2013l ) expand on this work, by adjusting their 
parameters to simultaneously match both the observed low- 
z Tf-band galaxy luminosity function and the HIMF. They 
find that this is not possible using standard models includ- 
ing the above assumptions, but it is possible if one does 
any of the following: 1) reduce Ecrit by about an order of 
magnitude, 2) allow the disk to maintain an exponential 
profile while stars are forming, i.e. allow gas from large radii 
beyond the star forming radius to lose angular momentum 
by some process and move inwards, 3) include some pro- 
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cess that preheats the gas before it can be accreted into 
a dark matter halo or have the formula by which one de- 
termines the expected gas accretion rate be different for 
small mass halos. The first option violates observations by 
a substantial margin, and furthermore for smaller, lower- 
met allicit^_gaJaxies_oiw_obs^^ to increase, not de- 
crease iBolatto. Wolfire. fc Lerovl2013l ). Our simulations as- 
sume a three dimensional critical star formation density of 
riH = 0.13 c m~^, which approximately matches the ob- 
served Tjcrit ISpringel fc Hernguistl 120031 ). and hence this 
is not likely to be the difi'erence between our models and 
the SAMs. The third possibility is also unlikely to be hap- 
pening in our models, since we do not explicitly add any 
sort of preheating, nor do we explicitly vary the dark mat- 
ter accretion rat es from the us ual ACDM expectatio ns (e.g. 
iNeistein fc Dekel l2008i : lFaucher^iguere et al.ll201lh . Hence 
if the arguments presented in ILu et al.l \2Qm apply to our 
simulations, then since we simultaneously match the HIMF 
and the stellar mass function, we suspect that angular mo- 
mentum transport processes, either physical or numerical, 
cause the gas in our simulated disks to move inwards. We 
intend to investigate this and other possibilities in future 
work. 

With the emergence of numerous major new radio facili- 
ties such as the Jansky Very Large Array (JVLA), MeerKAT 
(Karoo Array Telescope), and the Australian Square Kilo- 
meter Array Pathfinder (ASKAP), and eventually the SKA 
itself in the next decade, the future of H I studies looks 
promising. Simulations like the ones presented here, hope- 
fully improved substantially in the coming years, will be cru- 
cial for providing a context to these observations within the 
broader multi-wavelength landscape of galaxy evolution. At 
this point, it seems equally important both to probe the evo- 
lution of the H I content of galaxies and to go deeper around 
nearby galaxies to better understand how H I connects to the 
Cosmic Web. Such data promise to provide interesting new 
constraints on the key processes of galaxy evolution such as 
infiows, outfiows, and wind recycling. This paper represents 
a first step towards providing a comprehensive framework 
for interpreting the wealth of forthcoming H I data. 
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